module solve

    implicit none

contains

    pure elemental subroutine swap(a, b)
        real(kind=8), intent(inout) :: a
        real(kind=8), intent(inout) :: b
        real(kind=8) :: tmp
        tmp = a; a = b; b = tmp
    end subroutine swap

    !> code:blocks群友提供。
    !> 自动数组
    pure function solve1(a, b) result(x)
        ! n*n * n*m = n*n
        real(kind=8), intent(in) :: a(:, :)
        real(kind=8), intent(in) :: b(:, :)
        real(kind=8) :: x(size(b, 1), size(b, 2))

        real(kind=8) :: C(size(A, 1), size(A, 2) + size(b, 2))
        integer :: n, m
        real(kind=8) :: temp, temp2
        integer :: i, j, k, major

        n = size(a, 1)
        m = size(b, 2)
        C(:, :n) = a
        C(:, n + 1:) = b

        ! column major method
        do i = 1, n - 1
            ! 寻找第i列主元 find major position
            major = maxloc(abs(c(i:n, i)), dim=1) + i - 1
            !数据交换,cahce miss
            call swap(c(i, i:), c(major, i:))
            if (abs(c(i, i)) < epsilon(temp)) return
            temp = 1.d0/c(i, i)
            !计算
            do j = i + 1, n + m
                c(i, j) = c(i, j)*temp!归一化
                c(i + 1:n, j) = c(i + 1:n, j) - c(i + 1:n, i)*c(i, j)
            end do
            !第i列
            c(i, i) = 1.d0
            c(i + 1:n, i) = 0.d0
        end do

        temp = 1.d0/c(n, n)
        c(n, n) = 1.d0
        ! 回代
        do j = n + 1, n + m
            !归一化
            c(n, j) = c(n, j)*temp
            do k = 1, n - 1
                temp2 = c(n - k + 1, j)
                c(1:n - k, j) = c(1:n - k, j) - temp2*c(1:n - k, n - k + 1)
            end do
        end do

        x = c(:, n + 1:)

    end function solve1
    
    !> code:blocks群友提供。
    pure subroutine solve2(c, n, m)
        integer, intent(in) :: n, m
        real(kind=8), intent(inout) :: C(n, n+m)
        real(kind=8) :: temp, temp2
        integer :: i, j, k, major

        ! column major method
        do i = 1, n - 1
            ! 寻找第i列主元 find major position
            major = maxloc(abs(c(i:n, i)), dim=1) + i - 1
            !数据交换,cahce miss
            call swap(c(i, i:), c(major, i:))
            if (abs(c(i, i)) < epsilon(temp)) return
            temp = 1.d0/c(i, i)
            !计算
            do j = i + 1, n + m
                c(i, j) = c(i, j)*temp!归一化
                c(i + 1:n, j) = c(i + 1:n, j) - c(i + 1:n, i)*c(i, j)
            end do
            !第i列
            c(i, i) = 1.d0
            c(i + 1:n, i) = 0.d0
        end do

        temp = 1.d0/c(n, n)
        c(n, n) = 1.d0
        ! 回代
        do j = n + 1, n + m
            !归一化
            c(n, j) = c(n, j)*temp
            do k = 1, n - 1
                temp2 = c(n - k + 1, j)
                c(1:n - k, j) = c(1:n - k, j) - temp2*c(1:n - k, n - k + 1)
            end do
        end do

    end subroutine solve2

    !> walcs中求解线性方程组的例程，大致比欧拉的代码慢了4倍左右，且考虑了`pure`！
    !---------------------------------------------------------
    !用高斯消去法解方程组子程序,修改：参数h改为hh /by sun
    !---------------------------------------------------------
    pure subroutine Gauss(c, n, nm)
        !列选主元Gauss消去法解代数方程组
        !最后c(1:n,n+1:nm)中存放方程组的解
        implicit none
        integer(4), intent(in) :: n, nm        !n--增广矩阵的行数    nm--增广矩阵的列数
        real(8), dimension(n, nm), intent(inout) :: c        !增广矩阵
        integer(4) :: i, j, k, l, m, n1        !局部工作单元
        real(8) :: hh, g                       !局部工作单元
        n1 = n + 1
        do k = 1, n
            g = 0.0
            do i = k, n
                hh = c(i, k)
                if (abs(hh) <= abs(g)) cycle
                g = hh
                m = i
            end do                          !以上为选主元过程
            if (m /= k) then
                do j = k, nm
                    hh = c(k, j)
                    c(k, j) = c(m, j)
                    c(m, j) = hh
                end do
            end if                          !以上为行交换
            g = 1.0/g
            l = k + 1
            do j = l, nm
                hh = c(k, j)*g
                c(k, j) = hh
                if (l > n) cycle
                do i = l, n
                    c(i, j) = c(i, j) - c(i, k)*hh
                end do
            end do
        end do                                 !以上为消元
        do l = 1, n
            i = n1 - l
            do j = i + 1, n
                hh = c(i, j)
                do k = n1, nm
                    c(i, k) = c(i, k) - c(j, k)*hh
                end do
            end do
        end do                           !回代
    end subroutine Gauss
!--------------------------------------------------------

end module solve